*-------------------------------------------------------------------------------
* Heterogenous Relative Government Spending Multipliers in the Era Surrounding the Great Recession
* Forthcoming in the Review of Economics and Statistics
* by Marco Bernardini, Selien De Schryder and Gert Peersman (fall 2018)
*
* This file relates to figures 6(c) and 7(c)
*-------------------------------------------------------------------------------

////////////////////////////////////////////////////////////////////////////////
* 1) SETTINGS
////////////////////////////////////////////////////////////////////////////////
****************************
* housekeeping and loading *
****************************

* settings
clear all // clear memory 
cls // clear results window
set more off // do not display the "more" message
set matsize 300 // maximum number of variables in the model
set graphics off // do (not) show graphics in STATA
pause off // do (not) pause
adopath ++ PLUS

* globals
global expe "cnrecdeb_BP" // name of the experiment


* dataset
use "${path}panel_Q.dta", clear
keep stateid statename quarter gdprea govrea mor_dti // keep only analyzed variables
tsset stateid quarter


********************************************
* set parameters that govern specification *
********************************************

global p=4 //lag order
global hor=8 // maximum horizon
global band 68 90 // 68 (1SD), 90 (btw 1 and 2 SD), 95 (2SD)

global exerc_mar sdA sdB sdC sdD sdE sdF // exercises (margins)
global begsample=2005 // beginning of the sample (2005q1 earliest)
global endsample=2015 // end of the sample (2015q4 latest)

***************************
* construct raw variables *
***************************

* state-level NIPA variables
rename gdprea y			// state-level real GDP
rename govrea g			// state-level real government production
tsfilter hp cycle = y, smooth(1e6) trend(pot) // variable used to express changes in the same units
drop cycle

* economic states
gen lgdp=100*log(y) // log quarterly GDP
gen growth=(lgdp - L1.lgdp) // QoQ growth
gen rec=0
replace rec=1	if growth<0 & growth!=.
replace rec=0	if L.rec==0 & rec==1 & F.rec==0
by stateid: ipolate mor_dti quarter, gen(mor_dti_int)
gen debt = mor_dti_int // household debt-to-income

* fix analyzed sample
drop if  quarter<tq(${begsample}q1) // define the final sample (start)
drop if  quarter>tq(${endsample}q4) // define the final sample (end)

*****************************
* construct final variables *
*****************************

* LHS variables and reference for RHS lagged controls
foreach var in y g { // compute var_h for h=0,...,H
	forvalues j=0/$hor {
		gen `var'`j' = (F`j'.`var'-L.`var')/L.pot // DIFF(t+h - t-1) scaled by predet potential GDP (t-1)
	}
}

foreach var in y g { // cumulative LHS variables
	qui gen cum_`var'=0
	forvalues i=0/$hor {
		gen cum`i'_`var'=`var'`i'+cum_`var'
		replace cum_`var'=cum`i'_`var'
	}
	drop cum_`var'
}

gen sample=.
replace sample=1 if g${hor}!=. & y${hor}!=. & L${p}.y0!=. & L${p}.g0!=. // estimation sample (fixed along the horizon H)

* useful variables
by stateid: gen t = _n // count observations by ID
gen h=t-1 // horizon

* RHS variables (centered wrt deterministic variables -> only important for the nonlinear model)
** state variables
* define underlying state variables
gen nonlinA_v0pre=-(L1.lgdp-L5.lgdp)	// YoY negative growth (= 4-quarter backward MA of QoQ negative growth)
sum rec															if sample==1
local rec_freq=1-r(mean)
forvalues lambda=1(0.5)3 {
	local cc=`cc'+1
	gen nonlinA_v0_`cc'= (exp(`lambda'*nonlinA_v0pre))/(1+exp(`lambda'*nonlinA_v0pre)) if sample==1
	cumul nonlinA_v0_`cc', gen(nonlinA_v0_`cc'_ecd)
	label var nonlinA_v0_`cc'_ecd "{&lambda}=`lambda'"
	local gr `gr' scatter nonlinA_v0_`cc'_ecd nonlinA_v0_`cc' 	if sample==1 & nonlinA_v0_`cc'_ecd>0.5 & nonlinA_v0_`cc'>0.5, msize(tiny) ||
}
graph twoway `gr', xline(0.8,lcolor(red) lpatter(dash)) xline(0.9,lcolor(red) lpatter(dash)) yline(`rec_freq',lcolor(red) lpatter(dash)) ///
legend(rows(4)) xlabel(0.5(0.1)1,grid) ylabel(0.5(0.1)1,grid) name(calibration)
pause
gen nonlinA_v0= (exp(2.5*nonlinA_v0pre))/(1+exp(2.5*nonlinA_v0pre))
local staA bci
gen nonlinB_v0=(L1.debt)			// debt-to-income ratio
summarize nonlinB_v0											if sample==1, detail // 1=1SD
local nonlinB_v0grandmean = r(mean)
local staB deb
* center continuous state variables (relative variation)
qui reg nonlinA_v0 i.stateid i.quarter 							if sample==1
predict nonlinA_v1												if sample==1, r // centering
qui reg nonlinB_v0 i.stateid i.quarter							if sample==1
predict nonlinB_v1												if sample==1, r // centering
* z-scores continuous state variables (just to ease interpretation)
summarize nonlinB_v1											if sample==1, detail // 1=1SD
local nonlinB_v1stddev = r(sd)
replace nonlinB_v1=nonlinB_v1/`nonlinB_v1stddev' // standard idiosyncratic deviation in debt (i.e. total constant)
* define final states of the economy
gen sdA = 1														if sample==1 // neutral
gen sdB = nonlinA_v1											if sample==1 // additional effect of higher recession probability
gen sdC = nonlinB_v1											if sample==1 // additional effect of higher household debt
gen sdD = nonlinA_v1*nonlinB_v1									if sample==1 // interaction
local totstat sdA sdB sdC sdD // analyzed states
local ctrs `ctrs' sdB sdC sdD // include states as additional variables, except sdA to avoid multicollinearity with the constant term
pause
** lagged controls
foreach state in `totstat' {
	forvalues j=1/$p {
		capture drop l`j'g0 l`j'y0
		*
		qui reg L`j'.g0 i.stateid i.quarter if sample==1
		predict l`j'g0 if sample==1, r // centering
		gen l`j'g0_`state'=l`j'g0*`state'
		*
		qui reg L`j'.y0 i.stateid i.quarter if sample==1
		predict l`j'y0 if sample==1, r // centering
		gen l`j'y0_`state'=l`j'y0*`state'
		*
		local ctrs `ctrs' l`j'g0_`state' l`j'y0_`state'
	}
}

** shocks
qui reg g0 i.stateid i.quarter l1g0 l1y0 l2g0 l2y0 l3g0 l3y0 l4g0 l4y0 if sample==1 // fiscal rule + centering
predict shock if sample==1, r // residual
foreach state in `totstat' {
	gen shock_`state'=shock*`state' // state-dependent shocks							  
}
pause
* Pre-storage final objects
foreach ex_m in $exerc_mar {
** A) Cumulative multipliers
	*** 1-STEP (PEs and SEs)
	qui gen bmar`ex_m'=.
	qui gen semar`ex_m'=.
	foreach cint in $band {
		qui gen lomar`ex_m'_`cint'=.
		qui gen upmar`ex_m'_`cint'=.
	}

** B) IRFs
	foreach var in y g {
	*** 3-STEP (PEs and SEs)
		qui gen b`ex_m'_`var'=.
		qui gen se`ex_m'_`var'=.
		foreach cint in $band {
			qui gen lo`ex_m'_`cint'_`var'=.
			qui gen up`ex_m'_`cint'_`var'=.
		}
	}
}

gen Fsta = .

////////////////////////////////////////////////////////////////////////////////
* 2) DIRECT ESTIMATION
////////////////////////////////////////////////////////////////////////////////

forvalues i=0/$hor {
	capture drop cumfit_g_sdA cumfit_g_sdB cumfit_g_sdC cumfit_g_sdD
	local ord=`i'+1 // ord='i'+1 in ivreg2 is equivalent to ord='i' in xtscc (check "help ivreg2")
	global i=`i'

	** LP-IV
	qui reg cum`i'_g i.stateid i.quarter if sample==1
	predict cum`i'_gp if sample==1, r // centering
	gen cumfit_g_sdA = cum`i'_gp*sdA
	gen cumfit_g_sdB = cum`i'_gp*sdB
	gen cumfit_g_sdC = cum`i'_gp*sdC
	gen cumfit_g_sdD = cum`i'_gp*sdD
	
	ivreg2 cum`i'_y (cumfit_g_sdA cumfit_g_sdB cumfit_g_sdC cumfit_g_sdD = shock_sdA shock_sdB shock_sdC shock_sdD) `ctrs' i.stateid i.quarter if sample==1, dkraay(`ord') partial(`ctrs' i.stateid i.quarter)

	** storing
	local kpar = e(df_m)+1 // save # estimated parameters
	estadd scalar par = `kpar' // store # estimated parameters
	eststo	tabout_`i' // store regression output
	
	** exercises
	* base and marginal effects
	qui gen Fsta`i' = e(widstat)
	sum nonlinA_v0 if sample ==1
	local av_rec = r(mean)
	lincom cumfit_g_sdA - `av_rec'*cumfit_g_sdB
	qui gen bmarsdA`i' = r(estimate)
	qui gen semarsdA`i'= r(se)
	lincom cumfit_g_sdA + (1-`av_rec')*cumfit_g_sdB
	qui gen bmarsdB`i' = r(estimate)
	qui gen semarsdB`i'= r(se)
	lincom cumfit_g_sdC - `av_rec'*cumfit_g_sdD
	qui gen bmarsdC`i' = r(estimate)
	qui gen semarsdC`i'= r(se)
	lincom cumfit_g_sdC + (1-`av_rec')*cumfit_g_sdD
	qui gen bmarsdD`i' = r(estimate)
	qui gen semarsdD`i'= r(se)
	*
	lincom cumfit_g_sdB
	qui gen bmarsdE`i' = r(estimate)
	qui gen semarsdE`i' = r(se)
	lincom cumfit_g_sdD
	qui gen bmarsdF`i' = r(estimate)
	qui gen semarsdF`i' = r(se)
	
	qui replace Fsta = Fsta`i'	if h==`i' & Fsta`i'<80 & Fsta`i'!=.
	qui replace Fsta = 80		if h==`i' & Fsta`i'>=80
	foreach ex_m in $exerc_mar { // marginal effects
		qui replace bmar`ex_m'=bmar`ex_m'`i' if h==`i'
		qui replace semar`ex_m'=semar`ex_m'`i' if h==`i'
		foreach cint in $band {
			qui replace upmar`ex_m'_`cint'=bmar`ex_m'+invnormal((100+`cint')/2/100)*semar`ex_m' if h==`i'
			qui replace lomar`ex_m'_`cint'=bmar`ex_m'-invnormal((100+`cint')/2/100)*semar`ex_m' if h==`i'
		}
	}
	local tabout `tabout' tabout_`i'
}

////////////////////////////////////////////////////////////////////////////////
* 3) DECOMPOSITION
////////////////////////////////////////////////////////////////////////////////

forvalues i=0/$hor {
	foreach var in y g {
** 3-STEPS ESTIMATION
		ivreg2 `var'`i' shock_sdA shock_sdB shock_sdC shock_sdD `ctrs' i.stateid i.quarter if sample==1, dkraay(`ord') partial(`ctrs' i.stateid i.quarter)
	* individual responses
		lincom shock_sdA - `av_rec'*shock_sdB
		qui gen b`var'h`i'_sdA = r(estimate)
		qui gen se`var'h`i'_sdA = r(se)
		lincom shock_sdA + (1-`av_rec')*shock_sdB
		qui gen b`var'h`i'_sdB = r(estimate)
		qui gen se`var'h`i'_sdB = r(se)
		lincom shock_sdC - `av_rec'*shock_sdD
		qui gen b`var'h`i'_sdC = r(estimate)
		qui gen se`var'h`i'_sdC = r(se)
		lincom shock_sdC + (1-`av_rec')*shock_sdD
		qui gen b`var'h`i'_sdD = r(estimate)
		qui gen se`var'h`i'_sdD = r(se)
		*
		lincom shock_sdB
		qui gen b`var'h`i'_sdE = r(estimate)
		qui gen se`var'h`i'_sdE = r(se)
		lincom shock_sdD
		qui gen b`var'h`i'_sdF = r(estimate)
		qui gen se`var'h`i'_sdF = r(se)
		
	* impulse and cumulative responses
		foreach ex_m in $exerc_mar {
			qui replace b`ex_m'_`var' = b`var'h`i'_`ex_m'	if h==`i'
			qui replace se`ex_m'_`var' = se`var'h`i'_`ex_m' if h==`i'
			foreach cint in $band {
				qui replace up`ex_m'_`cint'_`var'=b`var'h`i'_`ex_m'+invnormal((100+`cint')/2/100)*se`var'h`i'_`ex_m' if h==`i'
				qui replace lo`ex_m'_`cint'_`var'=b`var'h`i'_`ex_m'-invnormal((100+`cint')/2/100)*se`var'h`i'_`ex_m' if h==`i'
			}
		}
	}
}

////////////////////////////////////////////////////////////////////////////////
* 4) OUTPUT
////////////////////////////////////////////////////////////////////////////////

** Table
esttab `tabout' using "${path}output\tab4b_${expe}.tex", ///
replace page(lscape) se nonum b(2) se(2) sfmt(2) obslast label nonotes ///
keep(cumfit_g_sdA cumfit_g_sdB cumfit_g_sdC cumfit_g_sdD) ///
nostar scalars("par Parameters")
local tabout

** Figure
gen zero=0

local marsdA "average effect" "{it:in expansions}" 
local marsdB "average effect" "{it:in recessions}"
local marsdE "average effect" "{it:difference}"
local marsdC "additional effect of debt" "{it:in expansions}"
local marsdD "additional effect of debt" "{it:in recessions}"
local marsdF "additional effect of debt" "{it:difference}"

foreach ex_m in $exerc_mar {
	* g
	twoway (rarea up`ex_m'_90_g lo`ex_m'_90_g  h if stateid==1 & t<=${hor}+1, ///
			fcolor(gs14) lcolor(gs13) lpattern(solid)) ///
		   (rarea up`ex_m'_68_g lo`ex_m'_68_g  h if stateid==1 & t<=${hor}+1, ///
			fcolor(gs10) lcolor(gs10) lpattern(solid)) ///
			(line b`ex_m'_g h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(orange)) ///
			(line b`ex_m'_y h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(none)) /// trick to get the same scale
			(line up`ex_m'_90_y h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(none)) /// trick to get the same scale
			(line lo`ex_m'_90_y h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(none)) /// trick to get the same scale
			(pcarrowi 0 0 0 8, lcolor(black) lpattern(shortdash) mcolor(black) msize(zero)), ///
			legend(off) name(fig_`ex_m'_g,replace) title("`mar`ex_m''",size(medium)) ///
			graphregion(color(white)) plotregion(color(white)) ///
			xla(0(1)${hor},grid labsize(medsmall)) /// yla(`la`ex_m'_g',grid labsize(medsmall)) ysc(`ra`ex_m'_g') ///
			ytitle("Size",size(medsmall)) xtitle("Horizon",size(medsmall))
	* y
	twoway (rarea up`ex_m'_90_y lo`ex_m'_90_y  h if stateid==1 & t<=${hor}+1, ///
			fcolor(gs14) lcolor(gs13) lpattern(solid)) ///
		   (rarea up`ex_m'_68_y lo`ex_m'_68_y  h if stateid==1 & t<=${hor}+1, ///
			fcolor(gs10) lcolor(gs10) lpattern(solid)) ///
			(line b`ex_m'_y h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(orange)) ///
			(line b`ex_m'_g h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(none)) /// trick to get the same scale
			(line up`ex_m'_90_g h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(none)) /// trick to get the same scale
			(line lo`ex_m'_90_g h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(none)) ///
			(pcarrowi 0 0 0 8, lcolor(black) lpattern(shortdash) mcolor(black) msize(zero)), ///
			legend(off) name(fig_`ex_m'_y,replace) title("`mar`ex_m''",size(medium)) ///
			graphregion(color(white)) plotregion(color(white)) ///
			xla(0(1)${hor},grid labsize(medsmall)) /// yla(`la`ex_m'_y',grid labsize(medsmall)) ysc(`ra`ex_m'_y') ///
			ytitle("Size",size(medsmall)) xtitle("Horizon",size(medsmall))
	* y/g
	twoway (rarea upmar`ex_m'_90 lomar`ex_m'_90  h if stateid==1 & t<=${hor}+1, ///
			fcolor(gs14) lcolor(gs13) lpattern(solid)) ///
		   (rarea upmar`ex_m'_68 lomar`ex_m'_68  h if stateid==1 & t<=${hor}+1, ///
			fcolor(gs10) lcolor(gs10) lpattern(solid)) ///
			(line bmar`ex_m' h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(midblue)) ///		
			(pcarrowi 0 0 0 8, lcolor(black) lpattern(shortdash) mcolor(black) msize(zero)), ///
			legend(off) name(figmar_`ex_m',replace) title("`mar`ex_m''",size(medlarge)) ///
			graphregion(color(white)) plotregion(color(white)) ///
			xla(0(1)${hor},grid labsize(medium)) /// yla(`la`ex_m'',grid labsize(medium)) ysc(`ra`ex_m'') ///
			ytitle("Size",size(medium)) xtitle("Horizon",size(medium))
}

* government spending
graph combine fig_sdA_g fig_sdB_g fig_sdE_g, c(3) name(fig_gA,replace) imargin(vsmall) ycommon ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero))
graph combine fig_gA, r(2) name(fig_g,replace) imargin(vsmall) ycommon title("government value added",size(medsmall)) ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero))

* gdp
graph combine fig_sdA_y fig_sdB_y fig_sdE_y, c(3) name(fig_yA,replace) imargin(vsmall) ycommon ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero))
graph combine fig_yA, r(2) name(fig_y,replace) imargin(vsmall) ycommon title("gross domestic product",size(medsmall))  ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero))

* IRFs
graph combine fig_g fig_y, c(2) name(fig,replace) imargin(vsmall) ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero)) xsize(7) ysize(9) ///
title("(c) continuous recession indicator",size(small))
graph export "${path}output\7_${expe}_gy.pdf", replace
graph export "${path}output\7_${expe}_gy.eps", replace
graph save "${path}output\7_${expe}_gy.gph", replace

* multipliers
graph combine figmar_sdA figmar_sdB figmar_sdE, c(3) name(fig_multA,replace) imargin(vsmall) ycommon ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero)) xsize(7) ysize(4.5)
graph combine figmar_sdC figmar_sdD figmar_sdF, c(3) name(fig_multB,replace) imargin(vsmall) ycommon ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero)) xsize(7) ysize(4.5)
graph combine fig_multA fig_multB, r(2) name(fig_mult,replace) imargin(vsmall) ycommon ///iscale(*.8) /// 
graphregion(color(white) margin(0 5 0 0)) plotregion(color(white) margin(zero)) xsize(7) ysize(4.5) ///
title("(c) continuous recession indicator",size(small))
graph export "${path}output\6_${expe}.pdf", replace
graph export "${path}output\6_${expe}.eps", replace
graph save "${path}output\6_${expe}.gph", replace

twoway		(line Fsta h if stateid==1 & t<=${hor}+1, lpattern(solid) lwidth(thick) color(midblue)), ///
			legend(off) name(figsta,replace) title("Fstat",size(medlarge)) ///
			graphregion(color(white)) plotregion(color(white)) ///
			xla(0(1)${hor},grid labsize(medium)) ytitle("Size",size(medium))  xtitle("Horizon",size(medium)) yla(0(20)80,grid) ysc(r(0 80))

graph export "${path}output\6_${expe}_Fstat.pdf", replace name(figsta)
graph export "${path}output\6_${expe}_Fstat.eps", replace
graph save "${path}output\6_${expe}_Fstat.gph", replace
